/*
 *  SCGlib : A C library for the spectral coarse graining of matrices
 *	as described in the paper: Shrinking Matrices while preserving their
 *	eigenpairs with Application to the Spectral Coarse Graining of Graphs.
 *	Preprint available at <http://people.epfl.ch/david.morton>
 *  
 *	Copyright (C) 2008 David Morton de Lachapelle <david.morton@a3.epfl.ch>
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 
 *  02110-1301 USA
 *
 *  DESCRIPTION
 *	-----------
 *    The kmeans_Lloyd function is adapted from the R-stats package.
 *	  It perfoms Lloyd's k-means clustering on a p x n data matrix
 *	  stored row-wise in a vector 'x'. 'cen' contains k initial centers.
 *	  The group label to which each object belongs is stored in 'cl'.
 *	  Labels are positive consecutive integers starting from 0.
 *	  See also Section 5.3.3 of the above reference.
 */

#include "igraph_memory.h"
 
#include "scg_headers.h"

int igraph_i_kmeans_Lloyd(const igraph_vector_t *x, int n, int p, 
			  igraph_vector_t *cen, int k, int *cl, int maxiter) {

  int iter, i, j, c, it, inew = 0;
  igraph_real_t best, dd, tmp;
  int updated;
  igraph_vector_int_t nc;
  
  IGRAPH_CHECK(igraph_vector_int_init(&nc, k));
  IGRAPH_FINALLY(igraph_vector_int_destroy, &nc);
  
  for (i = 0; i < n; i++) { cl[i] = -1; }
  for (iter = 0; iter < maxiter; iter++) {
    updated = 0;
    for (i = 0; i < n; i++) {
      /* find nearest centre for each point */
      best = IGRAPH_INFINITY;
      for (j = 0; j < k; j++) {
	dd = 0.0;
	for (c = 0; c < p; c++) {
	  tmp = VECTOR(*x)[i+n*c] - VECTOR(*cen)[j+k*c];
	  dd += tmp * tmp;
	}
	if (dd < best) {
	  best = dd;
	  inew = j+1;
	}
      }
      if (cl[i] != inew) {
	updated = 1;
	cl[i] = inew;
      }
    }
    if (!updated) { break; }

    /* update each centre */
    for (j = 0; j < k*p; j++) { VECTOR(*cen)[j] = 0.0; }
    for (j = 0; j < k; j++) { VECTOR(nc)[j] = 0; }
    for (i = 0; i < n; i++) {
      it = cl[i] - 1;
      VECTOR(nc)[it]++;
      for (c = 0; c < p; c++) { VECTOR(*cen)[it+c*k] += VECTOR(*x)[i+c*n]; }
    }
    for (j = 0; j < k*p; j++) { VECTOR(*cen)[j] /= VECTOR(nc)[j % k]; }
  }
  igraph_vector_int_destroy(&nc);
  IGRAPH_FINALLY_CLEAN(1);

  /* convervenge check */
  if (iter >= maxiter-1) { 
    IGRAPH_ERROR("Lloyd k-means did not converge", IGRAPH_FAILURE);
  }

  return 0;
}

